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ABSTRACT 



We use three-dimensional magnetohydrodynamical (MHD) simulations to 
study the formation of a corona above an initially weakly magnetized, isothermal 
accretion disk. The simulations are local in the plane of the disk, but extend 
up to 5 vertical scaleheights above and below it. We describe a modification to 
time-explicit numerical algorithms for MHD which enables us to evolve such 
highly stratified disks for many orbital times. We find that for initially toroidal 
fields, or poloidal fields with a vanishing mean, MHD turbulence driven by the 
magnetorotational instability (MRI) produces strong amplification of weak fields 
within two scale heights of the disk midplane in a few orbital times. Although 
the primary saturation mechanism of the MRI is local dissipation, about 25% of 
the magnetic energy generated by the MRI within two scale heights escapes due 
to buoyancy, producing a strongly magnetized corona above the disk. Most of 
the buoyantly rising magnetic energy is dissipated between 3 and 5 scale heights, 
suggesting the corona will also be hot. Strong shocks with Mach numbers ^ 2 
are continuously produced in the corona in response to mass motions deeper 
in the disk. Only a very weak mass outflow is produced through the outer 
boundary at 5 scale heights, although this is probably a reflection of our use of 
the local approximation in the plane of the disk. On long timescales the average 
vertical disk structure consists of a weakly magnetized (j3 ~ 50) turbulent core 
below two scale heights, and a strongly magnetized (f3 £ 1CT 1 ) corona which is 
stable to the MRI above. The largescale field structure in both the disk and 
the coronal regions is predominately toroidal. Equating the volume averaged 
heating rate to optically thin cooling curves, we estimate the temperature in the 
corona will be of order 10 4 K for protostellar disks, and 10 8 K for disks around 
neutron stars. The functional form of the stress with vertical height is best 
described as flat within ±2H Z , but proportional to the density above ±2H Z . 

For initially weak uniform vertical fields, we find the exponential growth of 
magnetic field via axisymmetric vertical modes of the MRI produces strongly 
buoyant sheets of magnetic energy which break the disk apart into horizontal 
channels. These channels rise several scale heights vertically before the onset of 
the Parker instability distorts the sheets and allows matter to flow back towards 
the midplane and reform a disk. Thereafter the entire disk is magnetically 
dominated and not well modeled by the local approximation. We suggest this 
evolution may be relevant to the dynamical processes which disrupt the inner 
regions of a disk when it interacts with a strongly magnetized central object. 



Subject headings: accretion disks - instabilities - MHD - turbulence 
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1. Introduction 

The observed spectrum of thin accretion disks in many systems suggests that the 
diffuse upper layers are heated to much higher temperatures than the midplane, i.e. that 
such disks have coronae. For example, disks around compact objects (Nandra & Pounds 
1994), CVs (Yi & Kenyon 1997), classical T Tauri stars (CTTS, Kwan 1997), and Herbig 
AeBe stars (Najita et. al. 1996) all show evidence for the existence of a corona. 

One simple explanation for the presence of a corona above a thin disk is that it is 
produced by external irradiation from the central object (or the central, hotter regions 
of the disk itself). A variety of authors have investigated the vertical structure of a thin 
disk atmosphere using radiative transfer calculations (Wade & Hubeny 1998; Dove et. al. 
1997; Haardt & Maraschi 1991) including the effect of irradiation from the central regions 
(D'Alessio et. al. 1998; Wood et. al. 1996). Such studies often find the diffuse upper layers 
of the disk are much hotter than the midplane. However, it should be noted that including 
a sophisticated treatment of the internal dynamics of the disk into such calculations is a 
formidable task, thus important ingredients such as the vertical profile of the heating rate 
due to accretion usually must be assumed using, e.g., the "a-disk" prescription (Shakura & 
Sunyaev 1973). 

A second possible explanation for coronae above accretion disks is that they are 
a direct consequence of the internal dynamics of the disk itself rather than a radiative 
transfer effect, much in the same way the solar corona is thought to be heated by dynamical 
processes lower in the atmosphere. This possibility was first suggested as a mechanism to 
explain the X-ray emission of Cygnus X-l (Liang & Price 1977, Galeev et. al. 1979) and 
of Seyfert AGN (Liang & Thompson 1979). In a detailed investigation by Stella & Rosner 
(1984), "magnetic buoyancy-driven" convection of strong flux tubes produced by shear in 
the disk was identified as the most promising source of the vertical energy flux needed 
to heat a corona. Although it was recognized that shear amplification of the field would 
result in stresses on the fluid that might give rise to angular momentum transport in the 
disk (Coroniti 1981; Sakimoto & Coroniti 1981), and that buoyancy might be the primary 
saturation mechanism for the field (Sakimoto & Coroniti 1989), modeling the coupled 
amplification, saturation, and angular momentum transport processes from first principles 
required multidimensional MHD models that were beyond the scope of these pioneering 
studies. Nonetheless, these early models were able to fit observations well enough to suggest 
that the internal dynamics of the disk might play an important role in determining the 
vertical structure. 

Strong magnetic fields (which have f3 < 1, where (3 = 8n P/B 2 with P the thermal 
pressure, and B the magnetic field strength) in a stratified accretion disk atmosphere are 
subject to both the Parker and magnetic interchange instabilities (e.g. Stella & Rosner 
1984). As an extension of the earlier work discussed above, the nonlinear evolution of 
these instabilities has been the subject of recent sophisticated numerical MHD modeling 
(Kamaya et. al. 1997; Basu et. al. 1997; Chou et. al. 1997; Matsumoto et. al. 1993; 
Kaisig et. al. 1990). Although such studies have helped clarify the effect of magnetic 
buoyancy in determining the vertical structure of the disk, they do not address the field 
amplification process itself, and therefore cannot investigate questions such as how strong 
flux tubes are formed in the disk in the first place. It is clear that field amplification, 
buoyancy, and angular momentum transport processes are all highly coupled, and therefore 
must be studied in concert for a complete picture. 
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In recent years, our understanding of angular momentum transport processes 
in accretion disks has progressed rapidly, in part due to the identification of the 
magnetorotational instability (MRI) in weakly magnetized accretion disks, and the 
realization of the fundamental role it plays in the process (see the review of Balbus 
& Hawley 1998 for a complete discussion). Time-dependent MHD simulations have 
demonstrated that the MRI saturates as MHD turbulence in three-dimensions (Hawley, 
Gammie, & Balbus 1995a, hereafter HGB1), and that this turbulence drives a non-helical 
dynamo which produces strong amplification of the magnetic field (Hawley, Gammie, & 
Balbus 1995b, hereafter HGB2; Brandenburg et al. 1995, hereafter BNST; Stone et al. 
1996, hereafter SHGB; Brandenburg 1998). Knowledge of the field amplification process 
allows us to address the problem of the formation of a disk corona from first principles. 
Tout & Pringle (1992) have presented an interesting analytic argument for the formation 
of a heated disk corona via a cyclic process involving the MRI, the Parker instability, and 
reconnection. However, the full details of such a (nonlinear) cycle require computational 
studies of the nonlinear evolution of the MRI in stratified disks. Current studies show 
that while buoyancy does result in a significant vertical flux of magnetic energy, it is not 
the dominant saturation mechanism for the instability: instead local dissipative processes 
dominate (BNST, SHGB). Still, these results hint at the possibility of the formation of a 
disk corona via buoyant rise of dynamo amplified magnetic field in that the scale height of 
the field was much larger than that of the gas, leading to a rapid decrease of (3 with height. 

To date numerical simulations of the MRI in stratified disks have been limited by 
the relatively small vertical extent that can be included in the calculations. This limit is 
imposed by the enormous range in dynamical timescales in a stratified disk: the density 
falls so rapidly with height in an isothermal disk (p oc exp(— z 2 )) that the Alfven speed 
becomes very large in the upper layers. For a time-explicit numerical algorithm (in which 
the timestep is limited in part by the crossing time of an Alfven wave across a grid zone), 
this makes calculations for many orbits of the disk impracticable with current computer 
resources. Thus, the simulations of BNST and SHGB were limited to two scale heights 
above the disk midplane, which was not enough to see the transition from a weakly 
magnetized turbulent core to a strongly magnetized corona which was stable to the MRI. 
Moreover, it is not clear what the effect of the boundary conditions applied at two scale 
heights may have on the field strength and geometry in this region. Ideally, one would like 
to apply outflow boundary conditions that allow mass, energy, and momentum to leave 
the computational domain, however, as discussed in SHGB, it is not possible to prescribe 
a self-consistent outflow boundary condition for a highly tangled, strong magnetic field 
without exerting anomalous stresses on the fluid near the boundary. Thus, in most of the 
simulations described by SHGB, periodic boundary conditions were applied at the vertical 
boundaries which, although they are clearly not a good representation of a real system, 
are not significantly less realistic than other viable choices, such as reflecting walls. A few 
simulations presented by SHGB which used different boundary conditions demonstrated 
that the vertical profile of horizontally averaged quantities near the midplane was not 
changed, indicating that the vertical structure is relatively independent of the boundary 
conditions. Moreover, the simulations of BNST used reflecting boundary conditions with 
a non-zero vertical component to the magnetic field (so that the turbulence in the disk 
generates a time-dependent current on the boundary); the resulting vertical structure of the 
disk was not significantly different than that reported by SHGB. 

In this paper, we describe an extension to time-explicit MHD algorithms which removes 
the timestep constraint imposed by the divergence of the Alfven speed in regions of very 
low density, and therefore allows numerical simulations over many orbital times of a disk 
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which spans many scale heights in the vertical extent. We use this extension to study 
the evolution of weakly magnetized, isothermal stratified disks over many orbital times 
using a grid which extends 5 scale heights above the disk midplane. We show that for a 
variety of initial field strengths and geometries, buoyantly rising magnetic field generated 
by MHD turbulence driven by the MRI near the midplane of the disk forms a strongly 
magnetized (/3 <SC 1) corona above two scale heights. Because it is strongly magnetized, the 
corona is stable to the MRI, and moreover the field becomes nearly force-free there (i.e. 
(V x B) x B 0). Thus, the field is no longer highly tangled near the upper boundary, 
and we are able to apply outflow boundary conditions appropriate to the system. These 
simulations allow us to investigate the formation and structure of a magnetized corona 
by calculating from first principles the coupled field amplification, buoyancy, and angular 
momentum transport processes. We are also able to measure properties such as the vertical 
profile of the stress and dissipation rate in the disk which may be useful in other contexts 
(such as detailed radiative transfer calculations of accretion disk atmospheres). Our models 
are still limited by the fact that they are local in the horizontal plane, which as we discuss 
may affect processes such as the generation of magnetocentrifugal winds from the disk 
surfaces. However, fully global models using the methods described here are underway. 

The paper is organized as follows: in §2 we describe our numerical methods, in §3 we 
present the evolution of a variety of weak magnetic field configurations in stratified disks 
which span many vertical scale heights, in §4 we discuss the implications of our results for 
the vertical temperature profiles in CTTS and compact objects, and in §5 we conclude. 



2. Numerical Methods 

We solve the equations of compressible MHD including the effects of resistivity in the 
"shearing box" approximation (see HGB1 for a detailed description). That is, we adopt 
a frame of reference which is local, Cartesian, and corotates with the disk at angular 
frequency Q corresponding to a fiducial radius R c . In this local frame, we use spatial 
coordinates x = (R — R ),y = {R<t> — &t), and z. Assuming a Keplerian rotation profile and 
expanding within a small volume about the fiducial radius (AR <C R ), the equations of 
motion can be expressed as: 

! + v.^o, (1) 

_Z + v • Vv = -VP/p + —(J x B) - 2Q x v + 3Q 2 xx - Q 2 zz, (2) 
at pc 

f = Vx(vxB) + ^VB, (3) 

where J is the current. The third and fourth terms on the RHS in the momentum equation 
[2] represent Coriolis and tidal forces in the rotating frame of reference, while gravitational 
acceleration in the vertical direction is represented by the last term. The final term in 
the induction equation [3] allows for field dissipation through an Ohmic resistivity with a 
constant magnitude 770- We assume isothermality and adopt an ideal gas equation of state 
with constant temperature, T, specified by our choice of the sound speed, c s (see below). 
By adopting an isothermal equation of state we assume that all thermal energy generated 
through the dissipation of kinetic and magnetic energies is immediately and efficiently 



- 6- 



radiated away and that there is an inflow of heat to adiabatically expanding regions (such 
as the buoyant magnetic flux tubes described in §3 below). Without a more sophisticated 
treatment of the thermodynamics of the gas, we cannot directly measure temperatures 
produced in our dynamical models (although by equating the energy dissipation rate 
measured in our simulations with an optically thin cooling rate, we can estimate equilibrium 
temperatures that might be reached in the corona, see §4). Moreover, if the vertical 
structure of the disk is significantly altered due to heating by the MRI, the rate of energy 
transport via buoyancy, and therefore the properties of the corona, may be affected. Also, 
we note that the assumption of isothermal evolution may increase the buoyancy of the flux 
tubes (see Torkelsson 1993) although this effect has not been studied in three dimensions. 
However, dynamical models which include a realistic treatment of radiation transport 
within an optically thick midplane and optically thin upper layers are beyond the current 
investigation. 

We note that the current (the second term on the RHS in equation [2]) has been 
modified to include Maxwell's displacement current (e.g. Jackson 1975); see the appendix 
for the definition of J as it is used here. Physically, the displacement current limits 
the speed of propagation of electromagnetic waves to the constant c. With the usual 
assumptions of ideal MHD, it can be dropped from the system. We keep it here to limit the 
propagation velocity of Alfven waves to a constant c = cu m , where we set cu m to be a value 
much less than the speed of light, but much larger than any relevant propagation speed 
within the disk. We further discuss the implications of including the displacement current 
in the evolution equations, and describe a technique for constructing a finite-difference 
approximation to the equation of motion [2] including it, in the appendix. 

We use an initial equilibrium state consisting of a density distribution which is constant 
on horizontal planes and follows a Gaussian profile in the vertical coordinate, 



where iff = 2c 2 s /Vt 2 is the thermal scale height. Our disk is thin; i.e., H z /R = 0.01. Initially 
all components of the velocity are zero except v y = Qx. The choice of Q is arbitrary in the 
shearing box formalism; for comparison with SHGB, we choose Q = 10~ 3 and 2c 2 s = 10~ 6 
(which yields H z = 1). Adopting p Q — 1, our initial pressure at the disk midplane is thus 
P Q = 5 x 10~ 7 . To seed the MRI, random perturbations with amplitudes of order 0.001c s 
are added to the velocities. 

We study the evolution of three different initial magnetic field geometries: 1) a purely 
toroidal field, 2) a Zero Net Z (ZNZ) field, and 3) a uniform B z field. The purely toroidal 
field is initialized only within \z\ < 2H Z such that j3 = 25 is constant within this region. 
The ZNZ field is generated by finite-differencing the following vector potential: 



p(x,y,z) = p Q e 



(4) 




A x — 



< y / 2P o //5(0) cos(27nr) cos(y) e -(l z ^ 2 ^) 4 : 2H Z < \z\ < AH Z 

: \z\> AH Z 



(5) 



where (3(0) = 25 is the plasma parameter at the disk midplane. This vector potential 
generates a field which has only a vertical component B z oc sin(27rx) cosy within two scale 
heights, and which is closed via loops between 2 and 4 scale heights. This configuration has 
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a null mean within —2 H z . Thus, within two scale heights, this geometry is equivalent to 
the ZNZ field studied extensively by SHGB. Finally, the uniform z field is a constant field 
initialized over the whole grid, with (3(0) = 25. 

Our simulations are computed using the 3D version of the ZEUS code (Stone & Norman 
1992a, 1992b). The Ohmic resistivity term is differenced using the method described in 
Fleming et. al. (1999). These authors found that Rcm = H z c s /r] < 10 4 ultimately quenches 
the MRI after many orbital times if the initial field configuration has a zero mean. In order 
to explicitly measure the reconnection and resistive heating rate in the corona, we have 
computed several models with Rcm = 20, 000, which, for the numerical resolutions used 
here, is large enough to resolve the resistive lengths but small enough not to completely 
damp the turbulence driven by the MRI. 

Our simulations consist of a grid of size 1H Z x 10H Z x 10H Z , with a range of —0.5H Z to 
0.5H Z in the x direction, to 10H Z in the y direction, and —5H Z to 5H Z in the z, direction. 
We have also computed a model in which the horizontal extent of the box is increased by 
a factor of two. Our standard resolution is 32 x 64 x 128 grid zones. We partition the 
zones in the vertical direction such that there are 64 zones in the region \z\ < 2H Z ; thus, 
inside this region our resolution is the same as SHGB. Above \z =2, the grid zones in 
the vertical direction increase logarithmically. We note that we have run one simulation 
with a grid which is nonuniform over the entire vertical domain to test that our choice of 
nonuniform zone spacing above ±2H Z does not affect the results reported here (particularly 
that it does not affect the vertical structure described below). There were no qualitative 
differences between the test simulation and our fiducial model (some quantitative differences 
did develop within ±2H Z due to differences in resolution). We have also run some higher 
resolution simulations to test the effect of resolution on our results. Such simulations are 
computationally intensive, prohibiting very high resolution runs. 

In the x and y directions, we adopt the periodic shearing sheet boundary conditions 
described in HGB1. In the z direction, we have implemented free (outflow) boundary 
conditions which allow mass, energy, momentum, and magnetic field to leave the grid. As 
described in Stone & Norman (1992a; 1992b), such outflow boundary conditions in the 
ZEUS code are implemented by using a constant (zero slope) projection of all the dynamical 
variables into the boundary zones. 



3. Results 

Table 1 lists the properties of the simulations discussed in this paper. Column 1 
gives the identifier for each run; column 2 gives the value of (3 at the disk midplane; 
column 3 indicates the value of Rtu\ and column 4 gives the value of q im used in the 
displacement current. Columns 5 and 6 list our numerical resolution and the ratio of the 
critical wavelength, A c (calculated using equation 15 of HGB1), at the midplane to the grid 
size, which is given by Ay for the toroidal field runs and by Az for the B z runs. When 
this ratio is greater than 5, the fastest growing wavelengths are well resolved (see HGB1). 
Column 7 indicates the total time of the simulation in orbits. Columns 8, 9, and 10 give the 
time and space averaged values of the magnetic energy and the azimuthal (y) components 
of the Maxwell and Reynolds stress for the regions \z\ < 2H Z (hereafter referred to as the 
"disk") and \z\ > 2H Z (hereafter referred to as the "corona"); these quantities have all been 
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normalized to the initial pressure at the midplane. If we adopt the Shakura-Sunyaev a 
convention, the time averaged value of a for each simulation is simply the sum of columns 
9 and 10. 



3.1. Toroidal Field Simulations 

We have performed toroidal field simulations with and without explicit resistivity (the 
runs with resistivity have been performed at two different resolutions); these are listed in 
Table 1 with the prefix "BY." Our fiducial run (BY1) is a standard resolution, nonresistive 
run with (3 = 25 in the magnetized region (\z\ < 2H Z ). This simulation has been run for 50 
orbital times. The critical wavelength for this run corresponds to 8.3lAy; thus initially it is 
well resolved. 

Figure 1 is a space-time plot showing the variation of horizontally averaged magnetic 
(top panel) and kinetic (bottom panel) energies as a function of vertical height and orbital 
time. The plot is constructed from the two dimensional function F(z,t) by averaging each 
quantity / in the horizontal plane: 

F(z,t) = 1 1 f{ ™'/' t)dxdy . (6) 
J J axay 

The plot of kinetic energy clearly shows the development of strong fluctuations associated 
with turbulence in the disk core after approximately 4-6 orbits driven by the nonlinear 
evolution of the MRI. We note that the growth rate and saturation time of the toroidal 
modes of the MRI observed here are in agreement with the toroidal field simulations of 
HGB1 (c.f. their figure 8) and SHGB (c.f. their figure 7). From the plot of the magnetic 
energy, periodic rise of buoyant magnetic field from the disk core over a few orbital times is 
clearly evident. The magnetic energy in buoyantly rising structures reaches a maximum at 
approximately 2.5 H z and then decreases. Thus, buoyancy appears to add both magnetic 
flux and heat (through the dissipation of magnetic energy) into the region z > 2.5H Z . Since 
kinetic energy is dominated by the mass distribution, comparison of the spacetime plots of 
the kinetic and magnetic energies indicate that the buoyantly rising field carries a small 
amount of mass into the coronal region. Only a small fraction of the mass and magnetic 
energy escapes the grid through the vertical boundaries; most remains in the corona. At 
the end of the simulation, the value of (5 at the midplane yields a ratio of A c /Ay = 6.14; 
thus, the fastest growing wavelengths are still resolved. 

To illustrate the vertical structure produced by the MRI, Figure 2 plots horizontally 
averaged values of the density, /3, kinetic and magnetic energies, and the azimuthal 
components of the Maxwell and Reynolds stress averaged over orbits 25 through 50, after 
the saturation of the MRI. We also list in Table 2 the volume-averaged values of each 
component of the magnetic and kinetic energies, the off-diagonal components of the Maxwell 
stress, the azimuthal component of Reynolds stress, and the density in the disk and coronal 
regions averaged over the same times. The graph of density shows that it is still roughly 
Gaussian. Flaring occurs in the wings due to the transport of mass into the corona via 
buoyant flux (from Figure 1 the expelling of fluid from the disk into the coronal region is 
periodic and correlated with the rise of magnetic flux). However, the inner regions of the 
disk remain almost unchanged; the total mass in the corona is only 1.03 % of the total mass 
within 2H Z . Table 2 confirms that the coronal densities are two orders of magnitude less 
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than in the disk midplane, and indicates that there are significant fluctuations in coronal 
region, with (Sp 2 ) 1 ^ 2 / (p) ~ 0.37. The plasma (3 is greater than one inside the disk and 
increases to nearly fifty at the midplane, roughly twice its original value. This means that 
the core of the disk remains unstable to the MRI throughout the simulation. Above ±2H Z , 
(3 is <^ 0.7, confirming that a magnetically dominated corona has indeed been produced. 

Above ±2.5H Z , (3 becomes small enough that A c exceeds 10H Z (the domain size in the 
azimuthal direction); this is why this region is stable to the MRI. Overall, (3 varies by more 
than three orders of magnitude from the corona to the disk midplane. The magnetic and 
kinetic energy plots also show that the corona is magnetically dominated. The magnetic 
energy in the coronal region exceeds the kinetic by more than an order of magnitude. Note 
that the magnetic energy in the corona has increased from zero to 0.042 P due to the rise 
of flux from the disk. It is approximately constant within the disk core. The magnetic 
energy peaks at roughly ±2H Z ; this same behavior was observed by SHGB (a further 
indication that the vertical magnetic energy profile described by these authors was not 
entirely a product of the periodic vertical boundary condition they used). The scale height 
of B 2 /(8tt) is roughly twice the density scale height. Referring to Table 2, we see that the 
azimuthal component of magnetic energy strongly dominates in both the disk and coronal 
regions. The radial component of the kinetic energy is about a factor of two greater than 
the vertical and azimuthal components in the disk. In the corona, however, pSv 2 /2P Q , which 
measures the energy in the velocity fluctuations, is nearly two orders of magnitude greater 
than the other components. 

In order to ensure the vertical profiles and volume averaged values shown in Figure 
2 and Table 2 respectively are not influenced by the time period over which they were 
averaged, we have reconstructed them with data which has been time-averaged over orbits 
15 through 50 (note that this period includes the strong transient peak in magnetic energy 
associated with the initial growth of the toroidal field MRI evident around orbit 10 in 
Figure 1) and over orbits 35 through 50 (beginning the averages well after the time when 
this initial peak has faded). We find that the time averages which begin at orbit 15 (see 
columns 8,9, & 10 of Table 1) and those which begin at orbit 35 lead to nearly the same 
vertical profiles and volume- averaged values as reported in Figure 2 and Table 2, indicating 
these values are not dominated by transients associated with any particular period of 
evolution. Moreover, in §3.1.2, we show that higher resolution simulations produce similar 
profiles, and in §3.2 we show that other field geometries also lead to similar profiles. This 
suggests the vertical structure reported here for our fiducial model is a general outcome of 
the evolution of the MRI in highly-stratified disks. 

From Figure 2, we see that the average of the azimuthal component of the Maxwell 
stress within ± 2 H z is a factor of 4.05 greater than the azimuthal component of the 
Reynolds stress, consistent with the toroidal field simulations of SHGB. This indicates 
that the turbulence and angular momentum transport within the disk are magnetically 
dominated. The profile of the Reynolds stress component is dominated by the density 
(although it is slightly flatter within the disk region), while that of the Maxwell stress 
shows a very different pattern. Note that the azimuthal component of the Maxwell stress 
is fairly constant within the disk region, rising only slightly (less than a factor of 2) at the 
disk midplane and edges (±2H Z ). This means that a is fairly constant in the disk; our 
simulations yield an average value of 0.027 for the a parameter within the disk core. Again, 
we note that we find approximately the same value of a (^ 0.03) whether or not we include 
the initial period when the MRI is saturating when we perform the temporal average. Thus, 
even though magnetic energy decreases slightly toward the end of the simulation (see figure 
5), significant angular momentum transport continues to take place. 



In the corona, we find that the azimuthal components of the Maxwell and Reynolds 
stresses drop sharply by 1-2 orders of magnitude over only 3 H z . The dramatic decrease 
in stress as we move into the coronal region indicates both that the transport rate has 
dropped and that the flow is much less turbulent here. The average value of the azimuthal 
component of the Maxwell stress in the coronal region is a factor of 5.0 smaller than in 
the disk (see Table 2), indicating the magnetic field is much less tangled. We find that the 
profile of the combined azimuthal stress as a function of vertical height is approximately 

flat within —2H Z and declines as e~ z above —2H Z . Therefore, the azimuthal component of 
the combined stress does not follow the density profile within the disk core. 

Figure 3 is a space-time plot of the off-diagonal components of the Maxwell stress 
tensor, as well as the azimuthal component of the Reynolds stress tensor. Note that all 
components except the B X B Z component of the Maxwell stress are largest in the region 
\z\ < 2H Z , confirming that the disk flow is turbulent while the coronal flow is smoother. 
However, although the magnitudes of the stresses are much lower in the corona than in 
the disk, both regions are characterized by large fluctuations {(5f 2 ) l ^ 2 /{f) ~ 1-0 (10.0 
for the B X B Z and the B y B z components) over the whole grid - see Table 2). The B X B Z 
component peaks between 2-3 H z , just above the core of the disk, in the region where 
the magnetic energy begins to dissipate. B z has a local maximum in this region as well, 
indicating that the field geometry is most strongly bent upwards by buoyancy between 
2-3 H z . The azimuthal component of the Reynolds stress confirms that the disk remains 
turbulent throughout the simulation. Diagonal streaks in the azimuthal components of 
the Maxwell and Reynolds stresses are clearly associated with the rising magnetic energy. 
Peaks in both (but particularly in the azimuthal component of the Maxwell stress) can be 
seen in the rising flux tubes. The B X B Z and B y B z terms are two (one) orders of magnitude 
smaller than the azimuthal component of the Maxwell (Reynolds) stress (see Table 2) and 
show a highly sporadic pattern. This suggests that pressure is a more dominant force than 
the off diagonal components of the stress in leading to disk flaring. The B y B z term is 
of particular interest as it is associated with the formation of a magnetocentrifugal wind. 
Our fiducial simulation produces an average value of 3.71 xl0~ 4 P o for B y B z above —2H Z , 
suggesting that only a very weak, unsteady wind is produced. We note, however, that 
the use of shearing sheet boundary conditions in the radial direction enforces a Keplerian 
rotation profile throughout the entire grid, even in the upper regions where the plasma is 
magnetically dominated. This may strongly affect the formation of a magnetorotational 

wind. Across the outer vertical boundary, average mass flux M = 2.21 x 10~ 4 OdC s , where 
aa is the surface density at the disk "surface," taken to be ±2H Z . 

Above ±2H Z , the azimuthal components of both Maxwell and Reynolds stresses drop 
sharply, decreasing by more than an order of magnitude. The low value of the Reynolds 
term in the coronal region is, however, more a reflection of the low coronal densities than an 
indication of smooth, quiescent flow (c.f. Table 2 which indicates that pSv 2 ~ 0.4 — 0.9P o in 

the corona). The coronal region is, in fact, in a highly dynamic state dominated by shocks. 
Figure 4 shows the vertical profile of the horizontally averaged rms velocity dispersion 
normalized by the sound speed, i.e., 



( 



(At;) 2 



1/2 



JJ(Sv x ) 2 + (5v y ) 2 + (Sv z ) 2 dxdy 
c 2 s J J dxdy 



(7) 




The profile shown in Figure 4 has been averaged over snapshots taken at seven equally 
spaced orbital times, between orbits 25 and 50. The core of the disk is dominated by 
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subsonic turbulence. The velocities become nearly supersonic above roughly 3 H z , the 
height at which the buoyant magnetic energy dissipates. In this region, shocks with Mach 
numbers as high as 2.9 are observed. These high velocities form as acoustic waves generated 
by the MRI carry part of the energy generated in the disk into the corona, steepening into 
shocks and dissipating there. The corona is not turbulent, but is highly dynamic. 

It is interesting to ask whether largescale magnetic field structure is generated in 
the corona as it becomes magnetized. Figure 5 is a plot of the relative strengths of each 
component of the field energy. With the onset of the MRI, the energy in all three field 
components increases. However, it is clear that within both the disk and the corona B 2 
dominates; it is more than an order of magnitude stronger than either of the other two 
components in both regimes. This is not surprising since differential rotation favors the 
production of B y . In the disk we find an ordering of By pa 15B 2 and B 2 pa 2B 2 ; this is 
consistent with previous three-dimensional studies of the MRI (e.g., HGB1, HGB2, and 
SHGB). Interestingly, the energy in the toroidal field is roughly equal in the disk and corona 
by the end of the run. This is due both to the efficiency of buoyancy in magnetizing the 
coronal region and to our shearing sheet boundary conditions which force Keplerian rotation 
and may thus be favoring the production of B y . After only 10 orbital times the corona 
is fully magnetized; its configuration is roughly that of a well ordered, toroidal magnetic 
field. In the corona, B 2 and B 2 Z are approximately equal. Differential vertical velocities 
therefore increase B z somewhat in the coronal region but are not enough to produce a 
strong B z field here; the toroidal component of the field is larger than the poloidal by more 
than an order of magnitude. This differs slightly from the results of SHGB, who report an 
ordering of B x pa 3B Z in both the corona and disk regions, suggesting the increase in the 
vertical domain of our simulations has allowed buoyant motions to develop more fully. We 
see largescale toroidal field configurations in both the corona and the disk, with small scale 
"wiggles" in the x and z directions. The slight secular decrease in the total disk magnetic 
energy is consistent with previous toroidal field studies; simulations on timescales longer 
than is possible here show the field never completely dies away. 



3.1.1. Dissipation of Energy 

Turning now to the resistive equivalent of our fiducial run (BYR1), we quantify the 
roles of dissipation and transport by examining the time rate of change of magnetic energy 
in the disk and coronal regimes. We use the fiducial resistive run for this analysis because 
we are able to explicitly measure resistive dissipation rates. In each of these volumes, 
magnetic energy is lost through a combination of the Poynting flux out of the volume 
and conversion into other forms of energy. Similarly, it is generated by conversion from 
other forms of energy (e.g. kinetic energy associated with the shear by the MRI) and the 
Poynting flux into the volume. If we let Sy denote the net generation of magnetic energy 
(embodying field generation by the MRI and energy loss due to conversion into other forms 
of energy), then we can write an equation of magnetic energy conservation as follows: 

S v = d fv(B 2 /8*)dV _ f {B 2 /An)v . + — <f (B • v)B • dS — [ V J 2 dV, (8) 
at Js An Js Jv 

where the last term on the right indicates losses due to resistive heating and the second 
and third terms represent the Poynting flux out of the volume. We can characterize the 



dissipation of energy throughout the volume by comparing the relative strengths of each of 
these components in the disk and coronal regimes. We will also compare with the heating 
rate due to compression and artificial viscosity, Rheat, given by: 




(9) 



where the components of q, given by 



qi = 4p(dvi/dxi) 2 , 



(10) 



are the artificial viscous pressure used in the ZEUS code to capture shocks (thus the viscous 
heating as we have defined it includes shock heating), and we will compare with the flux of 
kinetic energy, F kin : 



We will report each quantity in units normalized by the initial pressure, P Q , per orbit. 
Time averages are computed for orbits 15-50, after saturation of the MRI. In addition, we 
have plotted the instantaneous net magnetic energy generation rate (Sy) in Figure 6 as a 
function of orbital time in the disk (left panel) and coronal (right panel) regions. In each 
plot, Sy has been averaged spatially within ± 2 H z (left panel) and above ± 2 H z (right 
panel) and normalized by the initial pressure, P a . 

In the disk, Sy has a net positive time-averaged value of 0.10, the average magnetic 
energy generation rate is ~ —0.0005, and the average Poynting flux out of the volume 
is 0.10. Resistive heating is very small, with an average value of 0.0023. In the disk, 
the evolution of the Sy term (see Figure 6a) is dominated by largescale (SSy ^> (Sy)) 
fluctuations; these are related to the production of field via the MRI and to the loss of 
flux via buoyancy. The net positive time average value of Sy indicates there is a local 
increase in magnetic energy within the disk. This agrees well with the results of SHGB (see 
their figure 11); they also found the evolution of Sy to be dominated by fluctuations and 
report a net positive value within the disk core. (The major difference is that we do not see 
their initial spike in Sy; this is because we start with a toroidal magnetic field instead of 
a poloidal field). Our average value is more than an order of magnitude larger than that 
reported by SHGB, probably due to the fact that we measure Sy at two scale heights at 
the boundary between the turbulent core and stable corona, where buoyancy fluxes are a 
maximum (SHGB reported values measured at only one scale height). 

Comparing generation and loss rates in the disk, we find that energy generation via 
viscous heating dominates. Loss of magnetic energy due to resistive heating is quite low, 
representing only 0.71 % of the viscous heating rate. The average Poynting flux out of 
the disk represents 30.9 % of the viscous heating rate while the flux of kinetic energy is 
only 5.03 % and has a net negative average (which represents flow into, not out of, the 
disk). Thus, the rate of energy loss as flux tubes rise out of the disk is about one quarter 
of the total disk energy generation rate. This shows that the MRI regenerates field quickly 
enough to ensure a net increase in magnetic energy and to sustain turbulence and angular 
momentum transport for the duration of the run in spite of significant field loss due to 
buoyancy. We also note that at two scale heights, the amount of magnetic energy lost via 
buoyancy is more significant than the value reported by SHGB at one scale height. 

In the corona, the time average value for Sy is -0.012. The magnetic energy generation 
rate has an average of 0.0002, the net Poynting flux into the volume is 0.011 and the average 
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resistive heating is 0.00017. Comparing these values with those in the disk, we see that the 
situation in the corona is very different. The net Poynting flux rate is positive which means 
that more magnetic energy enters the corona than leaves it. However, the mean energy 
changes only slightly with time and the net value of the Sy term is negative. This indicates 
there is a local decrease in magnetic energy in the coronal region despite the net influx of 
field. We therefore conclude that the excess magnetic energy transported into the coronal 
region through the z = ±2H Z boundaries is dissipated there. Looking at Figure 6b, we see 
that Sy shows large fluctuations; these correlate very well with the rise of flux tubes into 
the volume. The correlation of the net flux into the corona with the fluctuations in the net 
magnetic energy generation rate (the Sy term) indicates that transport of energy dominates 
local processes in the corona. Therefore, although there is a net transport of MRI generated 
flux into the corona, this region remains stable: no magnetic energy generation occurs and 
the net magnetic energy transported into the corona is dissipated there. We note that 
SHGB also report a net loss of field in the coronal region, although, as before, their net 
value for Sy is an order of magnitude smaller. 

The energy in the net transported magnetic field is 60.1 % of the viscous heating in 
the coronal region; magnetic energy is thus the dominant form of heating here. However, it 
represents only 3.4 % of the viscous heating in the disk; therefore, it is only a small fraction 
of the overall energy budget. Resistive dissipation in the corona is surprisingly small. This 
is most likely due to the fact that the field structure is predominantly toroidal in both the 
disk and corona (as discussed above). Therefore, it is already fairly ordered as it rises into 
the corona. This also means that loss of magnetic energy in the coronal region is due more 
to its conversion into other forms of energy (e.g., kinetic energy) than to reconnection (as 
measured by r) ) . The flux of kinetic energy is comparable in magnitude to the Poynting 
flux, although its average value is negative, indicating a net flow both back into the disk 
region and out of the upper boundary of the simulation (roughly 10 % of the kinetic flux 
escapes through the upper boundary). Thus, energy enters the corona mostly in the form 
of magnetic energy. We conclude that the heating of the corona occurs via the conversion 
of magnetic energy into kinetic energy which is then dissipated by the viscosity. 



3.1.2. Effect of Resistivity and Resolution 

The evolution described above is typical of all toroidal field runs. The resistive runs 
show growth due to the MRI at slightly later orbital times, but otherwise follow the same 
qualitative evolution. As we increase the resistivity (i.e. from Rm = 20, 000 in run BYR1 
to Rm = 10,000 in run BYR2), the effect increases (i.e., growth occurs at increasingly 
later orbital times), but the general evolution does not change. The delay in growth with 
the addition of resistivity occurs because resistivity tends to diffuse gradients in the field, 
resulting in damping of the fastest growing wavelengths (Sano et. al 1998; Fleming et. al. 
1999). 

Run BYR3 is a high resolution version of the fiducial resistive run (BYR1); it was 
run on a grid containing twice the standard resolution. Again, the qualitative evolution 
remains the same as that described above. Comparing over the same time periods, the post 
saturation volume averaged values of magnetic energy increase by a factor of 1.3 in the disk 
region and 1.8 in the coronal region. We find the same ordering in the relative strengths 
of the individual components of magnetic energy in both regions regardless of resolution 
(see Figure 5). The azimuthal components of both the Reynolds stress and Maxwell stress 
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increase with resolution by factors of 1.6-1.5 (respectively) in the disk region. In the corona, 
the azimuthal component of the Reynolds stress remains nearly the same, increasing by 
only a factor of 1.2 while that of the Maxwell stress increases by a factor of 1.5. Thus, 
doubling the resolution has made important differences to some components of the energy 
and stresses, but in all cases this difference is less than a factor of two. The increases in 
magnetic energy and stress are probably due to an increase in the range of MRI unstable 
wavenumbers which are resolved in the simulation. Overall, the a parameter increases by 
a factor of 1.5 in the disk as we double the resolution. Most importantly, plotting the 
equivalent of Figure 2 for run BYR3, we find the profiles of all four quantities to be almost 
identical. (3 spans the same range in values from the disk midplane to the corona, and mean 
magnetic and kinetic energy values are the same in both regions as are mean values of the 
azimuthal components of the Maxwell and Reynolds stresses. The only noticeable difference 
in the graphs is in the magnetic energy and the azimuthal component of the Maxwell stress; 
these both dip to lower values at the disk midplane in the standard resolution runs instead 
of remaining approximately constant throughout the disk as they do in the high resolution 
run. Again, this indicates that we may not be resolving all of the unstable wavelengths 
in the disk midplane in the standard resolution simulations. Shock and resistive heating 
remain the same as in the standard resolution runs. Poynting flux also remains the same 
over the entire grid while the flux of kinetic energy increases with resolution by factors of 4.8 
in the disk and 2.0 in the corona. However, the velocity dispersion, which characterizes the 
turbulence of the flow, has the same profile and average values as in the standard resolution 
runs (see Figure 4). We find that the evolution of Sy in both the disk and the coronal 
regions is the same as that described above for the standard resolution runs. The average 
value of Sy in the disk increases by 15% in the high resolution run; the average coronal 
value remains unchanged. It appears that our standard resolution simulations capture the 
qualitative evolution very well, although they may underestimate some quantities, such as 
the a parameter. 

We have also varied the initial field setup to ensure that our choice of ± 2 H z as the 
initially magnetized region is not reflected in the results. We tried limiting this initially 
magnetized region to ± 1 H z as well as holding B constant within ± 2 H z , enforcing 
(5 =25 at the equatorial plane only; these are also listed in Table 1 as BY2 and BY3, 
respectively. None of these changes resulted in any significant differences with the general 
evolution described above. Finally, we repeated the fiducial run doubling the radial domain 
(run BY4). The small radial extent of our simulations may be limiting the off-diagonal 
components of the stress tensor and thus the production of magnetocentrifugal winds. Our 
small radial domain size may also be limiting transport rates and fluxes. However, doubling 
the radial domain size was not enough to produce any significant difference. Magnetic 
energy and magnetic and kinetic stresses actually decreased (by factors of 1.1-1.2) in both 
the disk and coronal regions when the radial domain was doubled. This leads to a slight 
decrease in the a value for the disk in run BY4 (by a factor of 1.2) compared to the value 
in the fiducial run. The Poynting flux and flux of kinetic energy remained the same at both 
the disk surface (± 2 H z ) and the outer boundary of the simulation box. There was also 
no change in the net rate of generation of magnetic energy (the Sy term) in the coronal 
and disk regions. We did not continue increasing the domain size due to the amount of 
computational time required. Clearly, the most fruitful means of investigating these issues 
is to perform fully global simulations which span a large range in radii. 
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3.2. Zero Net Z Simulations 

The Zero Net Z field simulations are identified in Table 1 by the prefix "ZN." Our 
fiducial simulation (ZN1) is a standard resolution, nonresistive run with /3(0)=25. The 
simulation was run for 50 orbital times. The critical wavelength for this run corresponds 
to 20.8 Az (within the disk midplane); it is well resolved. Note that in SHGB this field 
configuration was the fiducial model. In this paper, we adopt toroidal field runs as the 
fiducial model both because a toroidal field configuration is physically motivated by the 
differential rotation of the disk and because we have to close B z fieldlines in the coronal 
region, producing curved field lines there. (This was not an issue for the models presented 
in SHGB since they used periodic boundary conditions at ±2H Z .) The primary effect of the 
curvature of the field lines in the upper layers of the disk is to create tension forces which 
cause motions in the low density regions long before the MRI has saturated in the disk 
midplane. 

Despite the change in initial magnetic field configuration, the qualitative evolution of 
the ZNZ runs is very similar to that of the toroidal field simulations (described above). As 
expected (see HGB1, HGB2, SHGB), the MRI emerges earlier (after only 2-3 orbits) than 
in the toroidal field runs. Figure 7 is a space-time plot of run ZN1 showing the vertical 
variation of horizontally averaged magnetic energy (top panel) , the azimuthal component 
of Maxwell stress (2 nd panel), kinetic energy (3 rd panel), and the azimuthal component of 
Reynolds stress (bottom panel) versus orbital time. The MRI develops initially in the upper 
layers of the disk. By three orbits, however, the dominant modes are those at the disk 
midplane. As in the toroidal field runs, the MRI increases the magnetic energy within the 
disk region, and the magnetic energy rises buoyantly into the corona. Strongly magnetized 
regions rise with approximately the same slope (as measured in the space-time plots) as 
in the toroidal field runs and saturate at 3 H z , creating a magnetized corona. As in the 
toroidal field simulations, this process is periodic on a short (3-4 orbit) timescale. The 
major difference between run ZN1 and the toroidal field runs is that here the instability 
appears to die away at the midplane at the end of the run; the average value of f3 at the 
midplane at late times yields a ratio of A c /A,2 = 5.9, indicating that the critical wavelengths 
are just barely resolved. 

The main quantitative differences between the ZNZ and the toroidal field simulations 
lie in the magnetic quantities, which are uniformly smaller in the ZNZ runs. This is most 
likely due to the fact that the critical wavelengths are not well resolved in the midplane 
by the end of the ZNZ simulations. Figure 8 shows the horizontally averaged values of the 
density, /3, kinetic and magnetic energies, and the azimuthal components of the Maxwell 
and Reynolds stresses averaged over orbits 10-50, after saturation of the MRI. Comparing 
Figure 8 with Figure 2, we see that the range in density is approximately the same for 
both field configurations. The average coronal values of j3 are about the same, although 
/3 is 6 times larger in the midplane in the ZNZ run, so that it spans slightly more than 4 
decades over the vertical extent of the grid. The range in the azimuthal component of the 
Maxwell stress is about the same in both runs, although the absolute values are smaller 
in the ZNZ field runs. Probably the most notable feature of the profiles is a sharp dip at 
the midplane in both the magnetic energy and the azimuthal component of Maxwell stress. 
This is largely due to the strong transient burst of magnetic energy that occurs around 2 H z 
between orbits 30 and 40 (see Fig. 7). The kinetic energy follows the density profile closely 
and has about the same range as in the toroidal field runs. The azimuthal component of the 
Reynolds stress drops more sharply in the disk surface layers and corona compared to the 
toroidal field runs, although it spans approximately the same range. The disk a parameter 
for this field configuration is about 0.01, which is about a factor of 3 smaller than the 
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value found in the toroidal field runs. We assign no significance to the dip in the magnetic 
quantities at the midplane (since it is caused by a transient feature); neglecting the dip, the 
profile of the stress for these runs follows the form we found for the toroidal field runs; i.e., 
it is approximately flat within the disk core and is proportional to the density in the corona. 

In spite of its initially vertical structure, the field within the disk becomes predominantly 
toroidal due to differential rotation before beginning its buoyant rise and is thus toroidal 
in the coronal region as well (see Figure 9). Comparing Tables 2 and 3, we find the same 
relative ordering of the components of magnetic energy in both the disk and coronal regions 
as we did in the toroidal runs (see Figure 5 for comparison). The ordering of the kinetic 
energies in both the disk and coronal regions is also the same as in the toroidal field runs. 
As before, all quantities show large fluctuations. 

We find that energies and stresses are generally slightly lower in the ZNZ simulations 
than in the toroidal runs. Magnetic energies in the disk are factors of 3-4 lower in the 
ZNZ runs, and in the corona factors of 2-3 lower. Surprisingly, even the energy in the z 
component of the field are lower in the ZNZ field simulations, in spite of the initial vertical 
field configuration. Kinetic energies are lower by factors of 2-3 in the disk and 1-2 in the 
corona. Magnetic and kinetic stresses are also lower, with ZNZ values averaging 4 (2) 
times smaller in the disk for the magnetic (kinetic) quantities and 3 times smaller in the 
coronal region for both magnetic and kinetic quantities. Two notable exceptions are the 
—B X B Z term, which is nearly identical in the coronal region to the toroidal field value, and 
the coronal value of the —B y B z term, which is almost a factor of 10 smaller for the ZNZ 
field simulations. We find an a parameter of .0079 for our fiducial ZNZ run; this about 
a factor of 3 smaller than that cited for the toroidal field simulations (a = 0.027). The 
B y B z and B X B Z components of the stress tensor show a much stronger correlation with the 
rising magnetic flux tubes described above, indicating this field configuration may be more 
favorable to the production of ordered stress patterns which can fuel magnetocentrifugal 
winds; however, their time averaged values are still very low (7.2 xlO -5 P Q for the B y B z 
component) in this region. Magnetic fluxes are roughly a factor of 3 smaller than those 
measured in the toroidal field runs; kinetic fluxes are factors of 2-3 smaller. However the 
ratios of flux to energy for both magnetic and kinetic quantities are approximately the 
same; this may explain why the coronal field geometry at saturation is identical to that seen 
in the toroidal field runs. Overall, there are surprisingly few qualitative differences between 
simulations run with ZNZ and toroidal initial field configurations, although most quantities 
are slightly lower in the ZNZ runs. 

We note that we repeated run ZN1 with resistivity (run ZNR1), using Rm = 20K as in 
the fiducial resistive toroidal field simulation. We find that the effects of resistivity on the 
ZNZ simulations is the same as that discussed above for the toroidal field simulations. (See 
Fleming et. al. (1999) for a more complete discussion of the effect of resistivity on different 
initial field configurations.) 



3.3. Pure B z Simulations 

The evolution of a pure B z field geometry is interesting because it is relevant to a 
circumstellar disk threaded by a stellar magnetic field; i.e., a disk with a net vertical flux. 
Although far from the central star the stellar field is likely to be dipolar, in the local 
approximation adopted here we represent the stellar field as a uniform strength B z field. 
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These simulations are listed in Table 1 with the identifier "BZ." Our fiducial simulation 
(BZ1) is a standard resolution, nonresistive run with (3(0) = 25. This simulation has 
been run for 15 orbital times. However, we will only discuss the evolution up to 5 orbital 
times; beyond this point, the disk geometry is destroyed and the grid becomes magnetically 
dominated everywhere. The shearing sheet boundary conditions, which enforce Keplerian 
rotation everywhere, are then no longer appropriate for following the further evolution. 

Figure 10 is a space-time plot showing the vertical variation of horizontally averaged 
magnetic energy (top panel), the azimuthal component of Maxwell stress (2 nd panel), 
density (3 rd panel), and the azimuthal component of Reynolds stress (bottom panel) versus 
orbital time (along the x axis). It is immediately obvious from a comparison of Figure 10 
with Figures 1, 3, and 7 that the evolution of a uniform B z field is completely different 
from that of the other field geometries discussed in this paper. The difference is due to the 
development of the axisymmetric phase of the MRI known as the channel solution (Hawley 
& Balbus 1991; Goodman & Xu 1994). The channel solution develops at about 3 orbits; 
it is characterized by exponential growth of the magnetic field which causes the disk to 
separate vertically into horizontal planes (axisymmetric "channels" - see panel 3 of Figure 
10). Note that the magnetic energy and the azimuthal components of the Maxwell stress 
are large in the regions between the channels; these rise due to a combination of magnetic 
buoyancy and pressure. The uppermost of these channels rise out of the simulation box in 
less than 1 orbital time (these can just barely be seen rising from the low density surface 
layers of the disk). The rise of the channels which originate closer to the midplane is 
checked by the onset of the Parker instability, which causes undulations in the rising density 
channels. Mass slides along the bends in the channels back towards the midplane, reforming 
a semblance of a disk structure which oscillates about the midplane before reaching an 
equilibrium state. The simulation is magnetically dominated everywhere by the end of the 
run. 

Figure 11 is a snapshot of the vertical structure at 5 orbits, after the disk has been 
reformed. The quantities in Figure 11 have been averaged spatially in the horizontal 
directions. The plot of density shows that the bulk of the disk is at the midplane, although 
its initial symmetric structure has been perturbed. The plot of (3 indicates that both the 
disk and the coronal regions are magnetically dominated; /3 ^ 1.0 everywhere. We see from 
the plot of magnetic energy that the field has become very large at the midplane; it is of 
the same order of magnitude as the kinetic energy in the disk and is suprathermal over 
the entire grid. This is very different from the toroidal field runs (see Figure 2), where the 
disk field remained subthermal throughout the evolution. The strong magnetic pressure 
at the midplane drives the asymmetric density distribution. The plot of the azimuthal 
components of the stresses indicates that the Maxwell term dominates over the Reynolds 
term everywhere, and that the Maxwell term varies by at most a decade over the entire 
grid. Thus, after the disk reforms, matter is highly magnetically dominated in both the disk 
and coronal regions. The field structure is tangled and complicated, with By ~ -B^ ~ B 2 Z 
over the entire grid. 

We have also completed runs in which the initial density profile is a function of x 
(BZ2) and of x & z (BZ3). Such density distributions seed faster growing modes of the 
Parker instability, allowing it to curb the violent rise of the buoyant channels more quickly. 
In these simulations, mass loss due to the channels is less; however, the general evolution is 
qualitatively similar to that described above. 

We conclude that purely vertical fields quickly lead to a magnetically dominated disk 
everywhere. Simulations beginning with purely vertical fields were attempted by SHGB, 
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and similar evolution was noted, however because of the much smaller vertical domain 
used by these authors, the buoyant channels would hit the boundary at ±2H Z before the 
Parker instability could grow. Further simulations of this initial field geometry will require 
global models in order to capture the dynamics of the strongly magnetized regions properly 
(without enforcing Keplerian rotation). 



4. Discussion 

In this paper, we have shown that the evolution of a weakly magnetized stratified disk 
is characterized by the generation of magnetic field via the MRI and the loss of magnetic 
field through the buoyant rise of flux tubes (as long as the channel solution is not excited). 
The magnetic energy of buoyantly rising field saturates at about 2.5-3.0 H z and thereafter 
dissipates, serving to both magnetize and heat the corona (the region above ±2H Z ). The 
presence of a hot corona above an accretion disk has important observational consequences, 
as discussed in the introduction. It is therefore interesting to quantize the dissipation of 
total energy (E tot = E mag + E kin + E th ) as a function of height above the disk and to use 
this to calculate a thermal equilibrium temperature profile in the corona assuming the 
heating rate is balanced locally by optically thin cooling. In this section, we will estimate 
the coronal temperature distribution for two important cases: 1) a protostellar disk around 
a CTTS and 2) a disk around a neutron star, and compare our results with the coronal 
temperatures measured by observations. 

We first examine the case of a typical protostellar disk. In order to convert quantities 
measured in our simulations into physical units, we adopt the average CTTS parameters 
found in Beckwith et. al. (1993) (i.e., mass of disk O.O3M ; radius of disk k, 50AU). 
We assume our simulation box is located at R = 3R* (the inner edge of the disk) and that 
within two scale heights the disk is optically thick and has a temperature of 1000 K. The 
total dissipation rate as a function of vertical height, D, we approximate as: 

D(z) = R heat + ( I r]J 2 dxdy) x Az, (12) 
J s 

where Rheat is as defined in equation [9] (except here we integrate over a differential volume 
(corresponding to a horizontal plane with a height of one grid zone) instead of the entire 
simulation domain), and each quantity in the sum is averaged over horizontal planes at 
each vertical position in the numerical grid, and also averaged over time. We convert the 
total dissipation rate into an average per particle heating rate, T, in erg cm 2 s^ 1 via 

T = D/(eVn 2 ), (13) 

where e is a dimensionless heating efficiency and n is the average number density of particles 
within the volume V. We find that the shock heating term dominates the other dissipative 
components in our simulations; because it tends to be strongly localized, we set e = 10~ 3 . 
To be consistent with our definition of D, the volume V is the planar area Ax Ay times dz, 
where dz is the height of a vertical grid zone in units of H z . To compute the equilibrium 
temperature, we equate the per particle heating rate with an optically thin cooling rate, 
A, and solve for the temperature T. Using data from run BYR1, averaged over orbits 15 
through 50, we find the profile shown in Figure 12, where we have used the cooling rate 
given in figure 2 of Dalgarno & McCray (1972, hereafter DM) to yield the temperature 
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scale. Note that our assumption of a disk temperature at ±2H Z of T(±2H Z ) ~ 1000-K" is at 
least consistent with our result. In calculating the temperature at this point, we have used 
the curve in DM corresponding to the highest degree of ionization. This seems reasonable 
in the low density coronal regions; however, we note that different assumptions about the 
degree of ionization at 2H Z would lead to a different temperature for the same amount of 
dissipation. 

The thermal equilibrium temperature is of order 10 4 K in the coronal region. This is 
consistent with the warm temperatures Kwan et. al. (1997) need to account for the O I 
and Balmer lines of CTTS. For comparison, we apply our model to the case of 1548C27, 
using the stellar parameters given by Najita et. al. (1996). Asumming we are at their 
fiducial radius, we find a temperature ~ 9, 000-10, 000K for the inner disk atmospheric 
temperature. This is consistent with the temperature they derive to model the CO overtone 
emmision in this star. Since our simulations are local, we cannot meaningfully compute 
a radial temperature profile for the corona to compare to their data. We emphasize that 
the fluctuations in D about the mean value are large (factors of 4) which implies that 
the heating rate and peak temperatures in the corona will be highly time dependent. We 
find the temperature can vary by as much as two orders of magnitude on short (t <C t or b) 
timescales. Thus, our simulations show that a turbulent protostellar accretion disk produces 
a nonturbulent, magnetized, warm (T ~ 10 4 ) corona, with short timescale temperature 
fluctuations. 

We now turn to the case of an accretion disk surrounding a compact object. Adopting 
typical values of 1.0 M Q for the mass of the neutron star, 10 8 cm for a representative disk 
radius, and M — 1.6 x 10~ 9 M Q yr -1 as a typical accretion rate, we can use the formulae 
in Lovelace et. al. (1995) to derive a fiducial disk mass, velocity, and time. These allow 
us to convert the total dissipation rate into the per particle heating rate, T, as before. 
Using the neutron star parameters and comparing with the cooling rate given in figure 1 
of Raymond et. al. (1975), we find that a typical thermal equilibrium temperature in the 
corona is of order 10 8 K at 3H Z , which is in the X-ray range. Thus, we find the amount of 
energy generated by the MRI in the midplane of the disk and transported into the corona 
via buoyancy is sufficient to heat the corona to temperatures inferred from observations for 
disk coronae around compact objects. 

It is encouraging that the estimates presented in this section are consistent with 
temperatures suggested by current observational data as well as numerical calculations. 
However, we caution that there are many uncertainties in our estimates, most particularly 
in our choice of e, which, although motivated by the strong dominance of shock heating in 
our simulations, was admittedly ad hoc. If the heating were less localized, e could be as 
high as 1.0, which would change the temperature range significantly. 



5. Conclusions 

We have performed quasi-global MHD simulations of the dynamical evolution of 
vertically stratified accretion disks which demonstrate the generation of a magnetized 
corona. Initially weak magnetic fields in the core of the disk (below two scale heights) are 
amplified via MHD turbulence driven by the MRI. This field becomes buoyant and rises 
out of the disk. The magnetic energy density associated with buoyantly rising magnetic 
field peaks at about 2.5-3.0 H z and then dissipates, creating a magnetized, heated coronal 



-20- 



region above the disk core. We show that this process is highly time-dependent; throughout 
our simulations, the accretion disk continues to generate field via the MRI, which leads 
to the rise of flux out of the disk and its subsequent dissipation in the corona. On long 
time-averages, the vertical structure of the disk is a turbulent core which is unstable to 
the MRI surrounded by a highly dynamic, magnetically dominated ((3 ^ 10 _1 ), heated, 
stable corona. The largescale field configuration in both the disk and the coronal region is 
toroidal. Our simulations are the first to follow the highly coupled amplification, buoyant 
transport, and saturation processes which control the magnetic field, as well as the angular 
momentum transport process which controls the disk, from first principles. 

We have characterized the dissipation of energy as a function of vertical height for these 
simulations. We find that local processes dominate over transport in the disk region. This 
leads to heating of the disk by turbulence and to the continuation of angular momentum 
transport via the MRI throughout the simulations. In the coronal region, however, transport 
processes dominate because it is strongly magnetized (f3 < 1.0) and therefore stable to the 
MRI. Shocks with Mach numbers as high asMsi 2.9 are observed in the corona throughout 
the evolution; these are produced both by the steepening of MHD waves associated with the 
turbulence in the disk as they propagate into the low density corona, and via mass motions 
into the corona itself. By equating the time-averaged per particle heating rate in the corona 
with an optically thin cooling rate, we estimate the equilibrium temperature of the corona 
is of order 10 4 K for a typical accretion disk around a CTTS, and of order 10 8 K for a disk 
around a typical neutron star, although there is considerable uncertainty in these estimates. 
We note that the temperature can vary by as much as two orders of magnitude on short 
(t <C t or i>) timescales. The functional form of the stress as a function of vertical height is 
flat within ±2H Z but falls off proportional the density ( 1.6.. clS du Gaussian) above ±2H Z . 

We find that these results are general in the sense that they do not depend on 
resolution, resistivity (so long as the MRI is not quenched), nor initial field configuration 
(as long as the axisymmetric linear phase (the channel solution) of the MRI is not excited). 
They do depend on the initial field strength (if the field is strong, f3 > 1.0, and the disk 
becomes MRI stable). They may also be affected by the vertical structure of the disk. We 
assume an isothermal profile in all the simulations presented here; adiabatic disks may not 
be as strongly stratified and this may affect the buoyant transport rates. We are also limited 
by the use of the local approximation in the plane of the disk, which enforces Keplerian 
rotation everywhere, including in the strongly magnetized corona. This may suppress the 
generation of MHD winds. Extending our models to be global in all three directions is the 
natural next step in this study. 

We have also presented simulations in which the axisymmetric linear phase of the 
MRI is excited. Unlike the evolution described above, these simulations are dominated by 
the formation of density channels, which rise quickly out of the disk due to a combination 
of magnetic buoyancy and pressure. The rising density sheets are subject to the Parker 
instability, which results in mass sliding back to the midplane to reform a disk structure. 
This highly dynamic evolution may be relevant to the processes by which a stellar magnetic 
field can disrupt an accretion disk within the star-disk interaction region. 

We thank Steve Balbus for suggesting improvements to the Alfven speed limiter 
algorithm, and John Hawley and the referee (Ulf Torkelsson) for useful comments. This 
work was supported by the NSF through grant AST95-28299. 
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A. Large Alfven Speed Limiter via the Displacement Current 

An important aspect of this work is the extension of the simulation domain to include 
ten vertical scale heights (five above and five below the disk midplane). In time explicit 
codes such as ZEUS, the time step is constrained in part by the shortest crossing time of 
an Alfven wave across a grid zone. Thus, low density regions (such as those three to five 
H z above the disk midplane) where the Alfven speed becomes unmanageably large must be 
avoided. In the current simulations, we are able to work in the low density coronal region 
through the addition of the displacement current in the equation of motion [2], combined 
with a artificially low value for the speed of light c. This method was first used by Boris 
(1970) (see also the review by Brecht 1985). 

In practice, we set c = cu m (an arbitrary scaling speed) in the displacement current 
term in equation [2], which constrains VA/cu m to be ^ 1. We can thus choose cu m such 
that Ax/va is a reasonable value everywhere in the simulation domain. In this work, 
we choose cu m such that it is always approximately an order of magnitude greater than 
the (isothermal) sound speed; cu m = 5.66 xlCT 3 . This ensures that we allow a good 
dynamical range in the runs; i.e., the Alfven speed and the sound speed are allowed to differ 
significantly. It also makes our omission of relativistic terms introduced by the displacement 
current reasonable (i.e., VA/cn m ^ 1, but in general u/c& m <C 1, where v represents the 
other velocities in the simulation). We also note that we have tested the effect of this 
limiter by varying it over a range of values up to the value we report in this paper (compare 
columns 8, 9, & 10 of Table 1 for runs BY5, BY6, & BY1 and ZN2 & ZN1). As c Hm is 
increased, the only effect is an increase in the coronal kinetic energies and stresses. This 
is not unexpected, since we are limiting magnetic stresses, and therefore velocities, in the 
coronal region through the use of cu m . The amplitudes of the magnetic energy are not 
significantly changed. 

Our derivation proceeds as follows. Consider the momentum equation: 

<9v 1 
p- + pv-Vv=-VP + -(JxB), (Al) 
at c 

where vertical gravity and the shearing box rotational terms have been omitted. The 
displacement current enters through the definition of the current: 

J ^l VxB 4fl + ^ < A2 > 

where E is the electric field and £ + is the charge density. The last term represents the 
electric field force on the plasma. We can self-consistently include it in the definition of the 
current because it does not affect the energy balance in the code (e.g., v x B ~ E, and 
v • E = 0.0). The motivation for choosing this form for J will become apparent below. 

Assuming infinite conductivity^, E is given by the usual Ohm's law. Substituting in for 



2 We note that the terms introduced by including a finite conductivity are of order v/c <C 1; 
thus, we can self consistently assume infinite conductivity in Ohm's law. 
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J, we find: 

<9v 1 1 <9E 

P^T + pv • Vv = -VP + — (V x B) x B - — x B + — (v x B). (A3) 

ot An 4ttc ot c 

Substituting in for E and rearranging leads to: 

<9v B 2 1 1 8(v x B"l f + 

+ pv . w = -VP - V- + _(B ■ V)B + x B + |(v x B). (A4) 

We use the induction equation to bring B into the time derivative in the penultimate term 
on the right: 

p^+pv-W = -V/>-vf! + i(B.V)B + -^^^--^[IvQ 2 -(Q.V)Q] + ^Q, 

8t 87T 47T 47TC 2 <9t 47TC 2 2 C 

(A5) 

where Q = v x B. We wish to write the penultimate term on the right as the divergence of 
a tensor to guarantee conservation of angular momentum. Because V • Q ^ 0.0, we must 
add (and subtract) '^s-QV • Q' to do this. Then, using the equation of mass conservation 
to bring p into the time derivative on the left, we find: 

^-v^-v.H + ^M_ ?Q __i_ Qv . Q+ i: Q , (A6) 

where 

p a = ( p + + W " ( A7 ) 

P^ = -^[QiQi ~ (A8) 

Note that the extra term we added to the current (the last term in eq. (A6)) exactly cancels 

the non-conservation term introduced in writing V ■ Pd- This justifies the definition of the 
current in eq. (A2). 

We will now neglect the P D terms, which are of order v/c<^: 1. (Note that we can not 

simply drop the term - even though it is of the same order of magnitude as the Pd 
terms - because it is needed for conservation. This is especially important for the current 
study because of the key role angular momentum plays in the disk evolution which is driven 
by the MRI.) Canceling these terms, we then combine the time derivatives and write the 
mass density as a generalized tensor: 

B 2 . _ BjB>i , . 

to obtain 

^ = -V-P. (AID) 



and 
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The addition of the displacement current enters the momentum equation only via the mass 
density. We apply this change to the ZEUS code by replacing p with p* for the pressure 
(gas + magnetic) source terms. We then rederive the equations of MOC, following the 
method outlined in Hawley & Stone (1995; hereafter HS) for the 1.5 D system. There is no 
explicate change in the induction equation. For the momentum equation, if we ignore the 
pressure terms and the off-diagonal terms in the modified mass density, we are left with 



which simplifies to 



^ = -V-( P vv) + i-(B.V)B, (All) 



dp*(v x x + v y y) dpv x (v x x + v y y) t 1 D dB y 



dt dx 4-7T dx 



+ —B x —±y (A12) 



for the 1.5 D system. The x direction represents only the advection of momentum. For the 
y direction, after a little algebra and making use of the fact that terms of order v/c <C 1, 
we can write: 

n * dv v n *„ dv v , 1 R dB v / A 1 o\ 

p -W ~ ~ p Vx ~dx~ + 4^"a£" (A13) 

Comparing with equation (3.13) of HS, we see that the addition of the displacement current 
requires only that we replace p with p* = (1 + 47r B j> )p (where "c" has been replaced with 

" lim 

the limiter u cu m ") in the MOC equations. 

We note that some of the simulations reported in this paper do not include the £ + v 
term in eq. [A2]. We find that these runs show the same initial linear growth and saturation 
amplitude of the MRI as the simulations which do include this term. There is some 
indication that the long-term averaged values of the magnetic quantities (e.g., the magnetic 
energy) are slightly larger when this term is included. However, the large fluctuations which 
result from the chaotic nature of the MRI turbulence can lead to even larger differences in 
averaged quantities depending, e.g., on the time domain used in forming the average. Thus, 
averages must be performed over several hundreds of orbital times before secular differences 
in time averaged quantities can be considered meaningful. This suggests that the difference 
in the long-term average of the magnetic energy described above (where "long-term" refers 
to averages over < 50 orbits) is not significant. We stress that the presence of this term has 
no effect on the qualitative evolution discussed in the paper nor on the globally averaged 
quantities, such as the vertical profiles (see figs 2, 8, and 11). There are also no differences 
in the kinetic energies, stresses, or fluxes. The small increase in the time averaged value 
of the magnetic energy leads to a 3 % increase in the percent of energy dissipated in the 
coronal region (i.e., the Poynting flux increases as the magnetic field strength in the disk 
increases). However, we find that a different set of initial (random) perturbations can 
lead to differences of order 10 % in the solution (again due to the chaotic nature of the 
turbulence produced by the MRI) and therefore we do not assign any importance to this 
3 % difference. In spite of the lack of evidence that the inclusion of the £ + v term in eq. 
[A2] produces any real change in the results, we choose to include it in the fiducial toroidal 
(BY1) and ZNZ (ZN1) field simulations in order to strictly conserve angular momentum. 
We note that the other runs cited in this paper do not include this term. 
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Fig. 1. — Space-time plot showing the vertical distribution of the horizontally- averaged 
magnetic (top panel) and kinetic (bottom panel) energy as a function of time in the fiducial 
model BY1. The magnetic energy varies between 0.0 (black) and 0.846 P Q (red); the range 
in kinetic energy between 0.0 (black) and 1.24 P Q (red). 

Fig. 2. — Vertical variation of horizontally and temporally averaged quantities in model 
BY1. The time average is constructed between orbits 25 through 50. The horizontal average 
is taken over the f(f> plane as described in the text. 

Fig. 3. — Space-time plot of B x B y (top panel), pv x 5v y (2 nd panel), B x B z (3 rd panel), and 
B y B z (bottom panel) in model BY1. The maximum (minimum) value, colored red (purple), 
in each plot is 0.322 (-0.118) P Q , 0.132 (-0.055) P Q , 0.04 (-0.03) P Q , and 0.044 (-0.046) P Q , 
respectively. 

Fig. 4. — Horizontally and temporally averaged velocity dispersion in run BY1, normalized 
by the sound speed. 

Fig. 5. — Volume averaged By (solid line), B x (dashed line), B z (dot-dashed line), normalized 
by P Q , within the disk (left panel) and the corona (right panel) for model BY1. 

Fig. 6. — Net magnetic energy generation rate Sy in the disk (left panel) and the corona 
(right panel) in model BYR1. In both plots, Sy has been volume averaged and normalized 
byP . 

Fig. 7. — Space-time plot of magnetic energy (top panel), B x B y {2 nd panel), kinetic energy 
(3 rd panel), and pv x 5v y (bottom panel) in model ZN1. The magnetic energy varies between 
0.0 (black) and 0.443 P Q (red), the kinetic between 0.0 (black) and 1.175 (red), B x B y between 
-0.021 P a (purple) and 0.134 P Q (red), and pv x between -0.030 P Q (purple) and 0.078 P (red). 
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Fig. 8. — Vertical variation of horizontally and temporally averaged quantities in model ZN1. 
The time average is constructed between orbits 10 through 50. The horizontal average is 
taken over the f<j) plane as described in the text. 

Fig. 9. — Volume averaged By (solid line), B 2 X (dashed line), B\ (dot-dashed line), normalized 
by P , within the disk (left panel) and the corona (right panel) for model ZN1. 

Fig. 10. — Space-time plot of magnetic energy (top panel), B x B y {2 nd panel), density (3 rd 
panel), and pv x 5v y (bottom panel) in model BZ1. The magnetic energy varies between 0.51 
P Q (black) and 28.74 P Q (red), the range in density is 0.0 (black) to 8.69 (red), B x B y varies 
between -1.23 P Q (black) and 20.93 P Q (red), and the range in pv x 5v y is -5.1P D (blacK) to 
62.47 P (red). 

Fig. 11. — Horizontally averaged vertical structure at orbit 5 in model BZ1. The horizontal 
average is taken over the rcf> plane as described in the text. 

Fig. 12. — Vertical variation of coronal thermal equilibrium temperature for model BYR1 
for the case of a CTTS. 
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TABLE 1: 

Parameters of Numerical Simulations 



Run" 


0(0) 


Rm 




Grid 




A c /A 


Orbits 


/ B 2 v 6 
\ 8vrP(0) / 

(xlO- 2 ) 


/ -B X By I b 

\ 4ttP(0) / 

(xlO- 3 ) 


/ pV X SVy \ b 

\ P(0) / 

(xlO- 3 ) 


BY1 


25 


oo 


8 c s 


32 x 64 x 


128 


8.31 


50 


8.66/4.16 


23.6/4.33 


5.83/0.59 


BY2 C 


25 


oo 


8 c s 


32 x 64 x 


128 


8.31 


20 


6.34/1.80 


17.8/1.78 


4.36/0.394 


BY3 d 


25 


oo 


8 c s 


32 x 64 x 


128 


8.31 


20 


2.61/1.06 


7.12/1.43 


2.19/0.220 


BY4 e 


25 


oo 


8 c s 


64 x 64 x 


128 


8.31 


30 


5.87/3.16 


15.4/2.57 


4.37/0.389 


BY5 


25 


oo 


c s 


32 x 64 x 


128 


8.31 


50 


5.31/3.10 


13.7/3.72 


3.51/0.147 


BY6 


25 


oo 


4 c s 


32 x 64 x 


128 


8.31 


15 


9.41/6.72 


24.6/9.24 


6.03/0.629 


BYR1 


25 


20,000 


8 c s 


32 x 64 x 


128 


8.31 


50 


5.69/3.19 


14.3/3.17 


3.87/0.434 


BYR2 


25 


10,000 


8 c s 


32 x 64 x 


128 


8.31 


20 


5.55/2.98 


11.1/2.53 


2.85/0.340 


BYR3 


25 


20,000 


8 c s 


64 x 128 x 256 


16.61 


20 


10.9/6.17 


35.3/3.97 


8.58/0.582 


ZN1 


25 


oo 


8 c s 


32 x 64 x 


128 


20.8 


50 


2.20/1.28 


0.49/0.15 


0.23/0.021 


ZN2 


100 


oo 


c s 


32 x 64 x 


128 


10.38 


30 


0.372/0.313 


0.533/0.600 


0.269/0.045 


ZNR1 


25 


20,000 


8 c s 


32 x 64 x 


128 


20.8 


15 


0.486/0.131 


0.858/0.242 


1.29/0.046 


BZ1 


25 


oo 


8 c s 


32 x 64 x 


128 


20.8 


15 


179.2/101.3 


907.0/557.2 


56.9/12.0 


BZ2 


100 


oo 


C s 


32 x 64 x 


128 


10.4 


10 


88.1/42.6 


308.0/115.5 


35.5/6.99 


BZ3 


100 


oo 


C s 


32 x 64 x 


128 


10.4 


10 


103.4/50.1 


278.2/105.2 


31.4/2.60 



a Run labels: the first two letters give the initial field configuration (BY=toroidal field, 

ZN=zero-net-z field, BZ=pure z field). "R" identifies the resistive runs. 

b Volume and time averages begun at 15 orbits (BY6, where the average was begun at 10 

orbits and ZNR1, where the average was begun at 5 orbits). The first number represents 

the average within 2 H z ; the second represents the average for \z\ > 2H Z . For the BZ runs, 

these numbers are volume averages only, computed at orbit 5. 

c In this run, field was initialized only within —1H Z 

d In this run, (3 = 25 at the equator; B was held constant within —2H Z 

e In this run, the x direction was doubled without increasing the resolution: xe[— 1, 1]. 
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TABLE 2: 

Volume and Time Averaged Quantities in Run BYl 



Quantity f 


if) 


(5f 2 y/ 2a 


min f a 


max f a 




(xlO" 2 ) 


(xl0~ 2 ) 


(XlO" 2 ) 


(xl0~ 2 ) 




0.49/0.14 


0.23/0.09 


0.0/0.0 


1.0/0.58 


B 2 /8nP 


7.39/3.87 


2.15/2.58 


3.92/0.0 


12.6/8.93 


B 2 z /8nP 


0.26/0.18 


0.12/0.12 


0.0/0.0 


0.52/0.49 


-B x B y /AuP 


2.14/0.43 


0.97/0.27 


0.0/0.0 


3.80/1.47 


-B x B z /4nP 


0.01/-0.003 


0.10/0.04 


-0.09/-0.16 


0.16/0.20 


-B y B z /AnP 


-0.03/0.04 


0.05/0.05 


-0.30/-0.08 


0.21/0.15 


pv 2 x /2P 


1.11/0.12 


0.52/0.09 


0.0/0.0 


2.65/0.71 


pSv 2 /2P 


0.65/38.8 


0.28/20.0 


0.0/0.0 


1.22/88.3 


pv 2 /2P 


0.57/0.08 


2.43/0.06 


0.0/0.0 


1.23/0.29 


pVjV y /P 


0.53/0.06 


0.28/0.06 


0.0/-0.17 


1.36/0.31 


P 


64.0/0.71 


0.37/0.26 


63.6/0.314 


64.8/1.28 



a The first number represents the average quantity f within —2H Z ; the second number 

represents the average quantity f for \z\ > 2H Z . Averages have been taken over 25 orbital 
times, after saturation. 
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TABLE 3: 

Volume and Time Averaged Quantities in Run ZN1 



Quantity f 


(f) a 


(5f 2 y/ 2a 


min f a 


max f a 




(xl0~ 2 ) 


(xlO- 2 ) 


(xlO" 2 ) 


(xlO- 2 ) 




0.11/0.054 


0.059/0.036 


0.0/0.0 


0.23/0.19 


B 2 /8nP 


2.23/1.23 


1.06/0.67 


0.0/0.0 


3.89/1.96 


B 2 z /8nP 


0.060/0.052 


0.081/0.036 


0.01/0.006 


0.49/0.18 


-B x B y /4nP 
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a The first number represents the average quantity f within —2H Z ; the second number 

represents the average quantity f for \z\ > 2H Z . Averages were taken over 25 orbital times, 
after saturation. 
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